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ABSTRACT 

One of the parameters fitted by Doppler radial velocity measurements of extrasolar 
planetary systems is lv, the argument of pericenter of a given planet's orbit referenced 
to the plane of the sky. Curiously, the w's of the outer two planets orbiting Upsilon 
Andromedae are presently nearly identical: Auj = lo£, — ujc = 4?8 it 4?8(lcr). This 
observation is least surprising if planets C and D occupy orbits that are seen close to 
edge-on (sinic, siniD > 0.5) and whose mutual inclination 6 does not exceed 20°. In 
this case, planets C and D inhabit a secular resonance in which Aw librates about 0° with 
an amplitude of ~30° and a period of ~4 x 10^ yr. The resonant configuration spends 
about one-third of its time with |Acj| < 10°. If ^ 40°, either Auj circulates or the 
system is unstable. This instability is driven by the Kozai mechanism which couples the 
eccentricity of planet C to to drive the former quantity to values approaching unity. 
Our expectation that G < 20° suggests that planets C and D formed in a flattened, 
circumstellar disk, and may be tested by upcoming astrometric measurements with the 
FAME satellite. 

Subject headings: celestial mechanics — planetary systems — stars: individual (v An- 
dromedae) 

1. INTRODUCTION 

Upsilon Andromedae {v And) is a Sun-like star harboring at least three planetary compan- 
ions (Butler et al. 1999). The star's radial velocity variations are fitted approximately by the 
superposition of three Keplerian sinusoids, each of which takes the form 



V = K[cos{f + uj) + e cos uj] 



(1) 
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f is the true anomaly of a given planet, i is the unknown inclination between this planet's orbit 
plane and the plane of the sky, and the quantities G, M^, = l.SM©, m, a, and e take their usual 
meanings. Definitions and current fitted values of orbital parameters are listed in Table 1, as kindly 
supplied by G. Marcy and D. Fischer (2001). 

The quantity uj in equation (1) and Table 1 is a given planet's argument of pericenter referenced 
to the plane of the sky: 



u = sign[(n x e) • 1] arccos(n • e) G (— tt, tt] (3) 

where 



" — / • 

^1 - (S • i)2 

Let the origin be positioned at the barycenter of v And and take the sky plane to contain this 
origin. Then s is the unit vector pointing from the origin to the observer; I is the unit vector 
parallel to the planet's orbital angular momentum; n is the unit vector directed from the origin to 
the node of the orbit on the sky plane at which the planet is approaching the observer; and e is 
the unit vector pointing from the origin to the pericenter of the planet's orbit. Figure 1 illustrates 
the geometry. A single value of u corresponds to a volume of allowed parameter space swept out 
by vectors s, 1, and e. 

Curiously, the w's of the outer two planets C and D are presently nearly identical: Aa; = 
UJD — = 4?8 lb 4?8 (Icr). Let us define B = arccos (Ic • Id) to be the mutual inclination between 
the two orbit planes of planets C and D. If we assume for the moment that G = 0°, so that the 
orbits are co-planar and "co-rotating," then the observation that |Ati;| <C Irad implies the near 
perfect alignment of orbital pericenters. We may ask whether the alignment is merely coincidental, 
or whether a dynamical mechanism exists to lock the apsidal hues together. 

Precedents for apsidal locking exist in the Solar System. Eccentric planetary rings precess 
rigidly about Uranus and Saturn; the inner and outer elliptical edges of a given ring maintain the 
same line of apsides by a balance of forces due to the quadrupole field of the central planet, ring self- 
gravity, and interparticle collisions (Chiang k. Goldreich 2000). Of greater relevance to the case of 
V And is the observation that the periapsides of some asteroids librate (undergo small oscillations) 
about the apsidal line of Jupiter (Milani & Nobili 1984). Such asteroids inhabit a secular apsidal 
resonance in which their time-averaged apsidal precession rates match that of Jupiter. 
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Planets C and D oi v And might occupy a similar apsidal resonance. How might we use 
this possibihty to constrain the origin of these planets? Suppose that such a resonance operates 
if G = 0°, and suppose that it docs not operate otherwise. In the latter case, suppose that 
Auj circulates (cyclically runs the gamut between 0° and 360°), so that the probability of finding 
I Aw I ^ Irad at any given time is small. Given these assumptions, the observation today that 
|Aa;| <^ Irad might lead us to suspect that 9 = 0°. Co-planarity would argue, in turn, that 
planets C and D formed from a flattened, circumstellar disk. 

Even if the resonant dynamics are fully understood, this line of reasoning is suggestive at 
best. A quantitative, Bayesian assessment of the a posteriori probability distribution for Q given 
the observation that \Aio\ <C Irad requires knowledge of the a priori probability distribution for 
0. The latter distribution depends upon the unknown conditions of the formation of these two 
planets — the very object of our inferences. 

This paper examines the long-term evolution of the orbits of planets D, C, and to a limited 
extent B, within the multi-dimensional parameter space of present-day viewing and orbital geome- 
tries that are permitted by the Doppler observations. We seek to identify regions in this parameter 
space — in particular, all ranges of the mutual inclination 6 — that correspond to orbital evolutions 
for which |Aa;| <C Irad for all time; we suggest that the planets of v And reside in regions that 
are characterized by such apsidal resonances. The volume of parameter space to be surveyed is 
substantial; for instance, if 7^ 0°, then a;/? = loq does not even necessarily imply that ec, e/j, Ic, 
and Id lie in the same plane. Equality of w's merely reflects the equality of angles which may well 
be referenced to different nodes and which may therefore have no direct physical relation. 

Previous explorations of this parameter space by Rivera &: Lissauer (2000ab) and Stepinski, 
Malhotra, Sz Black (2000) have indeed uncovered the existence of a secular apsidal resonance 
between planets C and D in the = 0° case. Furthermore, considerations of dynamical stability 
presented by these authors limit values of sinz in the = 0° case to be greater than ~0.3. Stability 
considerations also argue against > 60° (Stepinski et al. 2000). 

Our analysis extends and improves upon these previous results in a number of respects. First, 
we take explicit account of the dependence of to on viewing geometry. To our knowledge, the only 
study to mention this dependence is that of Rivera & Lissauer (2000a), who incorporate it in their 
calculations of co-planar, counter-rotating orbits (0 = 180°). Second, we systematically explore 
all of the parameter space spanned by possible orbital configurations of planets C and D and by 
possible lines-of-sight; previous studies explore only a fraction of this space. Third, we focus on the 
evolution of Au as a means of distinguishing between likely and unlikely orbital configurations. As 
we shall see, the existence of a secular apsidal resonance characterized by small libration amplitude 
in select regions of parameter space suggests tighter constraints on than considerations of stability 
alone. Finally, the Doppler-fitted values of orbital parameters that we employ are improved over 
those used in previous studies due to the increased number of data points {N = 189 radial velocity 
measurements as of January 2001; Marcy & Fischer 2001). 
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In §2, we take G to be zero and study the secular evolution of the orbits of planets C and D while 
neglecting the presence of the innermost planet B. We recover analytically the secular resonance 
found numerically by Rivera & Lissauer (2000ab) and Stepinski et al. (2000) and describe its 
physical character. In §3, we allow 7^ 0° while continuing to ignore planet B. All viewing/orbital 
configurations for which Ao; potentially librates about 0° are identified. In §4, we ask which of 
the resonant configurations found in §3 are likely to remain resonant with the inclusion of planet 
B. Section 5 marshals our results to argue that the planets of v And most likely originated in a 
circumstellar disk. 



2. Planets C and D: 6 = 0° 

In the co-planar, co-rotating case, the individual cj's represent longitudes of pcricentcrs ref- 
erenced to a common nodal vector fi and a common orbital plane. We ask in this case whether 
Aw circulates or librates, and investigate how the answer depends on uncertainties in the orbital 
parameters, including sinz. 



2.1. Analytic Description 

Long-term variations in the eccentricities and apsidal longitudes of the outer two planets are 
qualitatively well described by the classical secular solution of Laplace-Lagrange (L-L; see, e.g., 
Murray & Dermott 1999). Let q = m^o/mc = 2.03±0.05 (la), a = ac/ao = 0.3238±0.0009 (Icr), 
and let (3 = h^^^^{a) /h^^^^{a) = 0.476/1.19 = 0.399 ± 0.001 (la) be the ratio of Laplace coefficients 
(Brouwcr & Clcmcncc 1961). These quantities are evaluated using the orbital parameters of Marcy 
& Fischer (2001). Then 



ecexpicjc = ec-t-expz(5+t-|-7+) -|- ec- exp i(5_t -|- 7_) (5) 
CDexpiwo = eD+expi(c/+t-|-7+) -I- eD-expi(5_t-)-7_) (6) 

where the frequencies 5+ and g-, and ratios eD+/e.c+ and e£)_/ec-, are constants specified by the 
masses and secularly invariant semi-major axes of the two planets: 



(8) 



/ CD \ _ q - =F V(g - V^Y + 4gi/a^ 
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Here Pc is the orbital period of planet C. These equations neglect terms of order and {m/M^)^. 
The Laplacc-Lagrange description is valid for small values of the mutual inclination Q provided u; 
is replaced by u), the longitude of pericenter referenced to the invariable plane (the usual "dog-leg" 
angle defined in, e.g., Murray & Dermott 1999). 

In writing equations (5)-(8), we have omitted terms due to the third, innermost planet B. 
Provided planet B remains on an orbit that is more than a factor of 10 smaller than the orbits 
of the outer two bodies, its time-averaged potential acts mainly as a static quadrupole moment. 
This quadrupole could add corrections up to order {ruB / 'mc){aB / acY ~ ^ few percent to our 
expressions for the eigenfrequencies (7) and eigenvector amplitude ratios (8). In this section and 
§3, we neglect the effects of planet B for clarity of presentation and computational ease. We 
incorporate quantitatively the effects of planet B in §4. 

The remaining four constants ec+/ec-^ a-^d e/)_|-/e£i_ require specification of the ini- 
tial eccentricities and apsidal longitudes. If initial conditions are such that \ec+/c-C-\ ^ 1 <^iid 
\eD+/eD-\ ^ 1, then equations (5)-(6) imply that the two planets precess in lockstep at frequency 
g- with fixed eccentricities and aligned pericenters. The present-day orbital parameters of v And 
are such that most, but not all, of the power is in the "— " eigenmode. Setting uc = 0°, ujd = 4?8, 
ec = 0.25, and eo = 0.34 at t = 0, we obtain with a numerical root-finder: 



The evolution described by equations (5)-(9) is plotted in Figures 2abc. On average, both too and 
ioc advance at the same rate In addition, however, the admixture of power from the mode 
{ec+/ec- is not substantially smaller than unity) implies that ujc librates about at a faster 
frequency The argument of pericenter of planet D librates less since that planet carries the 

bulk (~75%) of the angular momentum in the system. 

The apsidal libration amplitude in L-L theory is a function of the initial individual eccentricities 
(eco and e^o)) the initial difference in apsidal longitudes (Aa;o), a, and q. For planets C and D, most 
of the uncertainty in the L-L libration amplitude arises from uncertainties in the fitted eccentricities. 
Figure 3a exhibits the variation of L-L libration amplitude with eco and euo, for fixed a = 0.324, 
Kd/Kc = 1.19, and Aujq = 0°. Allowance is made for the fact that for fixed a and Ku/Kc, 



the observed q = mD/mc oc y 1 — e^o/y ~ ^co according to equation (2). We take Awq to 
be 0° because this approximation permits the immediate solution of equations (5)-(8). Variations 
in libration amplitude with Auq G [— 4?8, 14°] (2cr confidence interval) are typically less than 10% 
(other parameters being held fixed); these variations are small compared to the ~50% variation 
in libration amplitude associated with the 2a confidence region in fitted eccentricities. Variations 
in libration amplitude associated with uncertainties in a and Ku/Kc arc negligibly small, at the 
level of ~0.3% and ~3%, respectively. Since the L-L libration amplitude is sensitive to the mass 
ratio but not to the individual masses, it is independent of sini in the co-planar case. 



7+ = -10°, 



7_ = 4?3, eD+/eD- = -0.030, ec+/ec- = 0.41 . 



(9) 
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(a) L-L Libration Amplitude (deg) 




0.209 0.224 0.239 0.254 0.269 0.284 0.299 

Eccentricity C 



(b) Max Acj (deg) / = 0° / sin i = 1.0 




0.209 0.224 0.239 0.254 0.269 0.284 0.299 

Eccentricity C 
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The L-L amplitude associated with the best-fit eccentricities is 25°. The corresponding fraction 
of time that the system spends with \Alo\ < 10° equals Pume = 0.26. Within the 2a error ellipse 
of fitted eccentricities, the amplitude ranges from 14° {Pume = 0.52) to 38° {Pume = 0.16 

Wc conclude from this analytic study that if planets C and D occupy co-planar, co-rotating 
orbits, then they inhabit a secular resonance for which Au librates about 0° with an amplitude of 
25° ± 6° (la). In other words, if 6 = 0°, the observation today that |Aa;| < 10° should not be 
regarded as surprising, since the best-fit system spends about one-quarter of its time in this state. 
These libration parameters change slightly as we incorporate more realistic effects in subsequent 
sections (e.g., the presence of planet B is included in §4). 

2.2. Numerical Solution 

We verify the analytical considerations above by numerically integrating the orbits of planets 
C and D using the SWIFT software package developed by Lcvison &; Duncan (1994). We employ 
their mixed variable symplectic integrator, which is derived from the algorithm invented by Wisdom 
& Holman (1991). As before, the effects of planet B are neglected. The timestep is sufficiently 
small {At = 3 days) that the total energy in the system is conserved to within 1 part in 10^ after 
t = 10^ yr. Initial conditions are taken from the best-fit parameter values in Table 1. The starting 
epoch of our integrations was chosen arbitrarily to be JD 2450160.1 days. Since the effects in 
which we are primarily interested arc secular (orbit-averaged), our results should be qualitatively 
insensitive to choice of starting epoch. We have verified that none of our conclusions changes by 
choosing an alternative starting epoch of JD 2451160.1 days. Moreover, since substantial changes 
in the osculating orbital elements of planets C and D occur only over secular timescales that are 
long compared to the duration of the radial velocity observations, the errors in the fitted orbital 
parameters in Table 1 introduced by the assumption of fixed Kcplerian ellipses are no more than 
several percent (Laughlin 2001; see also Laughlin & Chambers 2001). 

As seen in Figures 2def, the behavior of Au with time for sin i = 1 computed using SWIFT is 
similar to that predicted by L-L. For smaller values of sini G [1;0.3] and correspondingly greater 
planetary masses, the system remains in the secular apsidal resonance, but short period terms in 
the disturbing potentials effectively increase the amplitude of libration. Thus, the fraction of time 
for which |Ati;| < 10° decreases from Pume = 0.30 to 0.20 as sin? decreases from 1 to 0.3. For 
sini < 0.3, we find, as do Rivera & Lissauer (2000b) and Stepinski et al. (2000), that the system is 
unstable. 

Figure 3b plots the maximum value of Aa; over 10^ yr as a function of the initial eccentricities 
for sini = 1. The contours are similar to those predicted by L-L, with deviations introduced by 
short period terms and higher order secular terms. The range in Pume within the 2a error ellipse 
of fitted eccentricities is [0.20,0.54]. 
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3. Planets C and D: 9 7^ 0° 

How does the secular evolution of Ao; change with mutual orbital inclination, O? In particular, 
for what ranges of can planets C and D occupy an apsidal resonance? We answer these questions 
using Monte-Carlo sampling techniques and numerical orbit integrations. As in §2, we neglect the 
effects of planet B. We begin by delineating the parameter space spanned by viewing geometries and 
orbital configurations of planets C and D that are allowed by the Doppler observations. Subseqiient 
subsections describe (1) the method by which we sample the allowed volume of parameter space 
and (2) the dynamics that play out within this volume. 

3.1. Geometrical Considerations 

Let the set of "initial orientations" refer to the set of {Ic, Id, ec, e^i, s} such that eu • Id = 
ec ■ Ic = and Id ■ Ic = cos Go- The subscript "0" denotes the initial or present-day value; as 
we shall see, the mutual inclination varies with time in most circumstances. These orbital and 
viewing orientation vectors sweep out a volume in the following five dimensions. Without loss of 
generality, we erect a set of Cartesian coordinate axes x = ec, z = Ic, and y = z x x. Then the 
five dimensions of initial orientation space are Qq,^Iq,xq — the initial inclination, initial longitude 
of ascending node, and initial argument of pericenter, respectively, of the orbit of planet D upon 
that of planet C — and 0, S — the azimuthal and polar angles of s. In the absence of observational 
constraints, the angles f^OjXOj^ € [0, 27r), while Qo,S G [0,7r]. 

The observed present-day values of ujc and wd carve out a subset in initial orientation space 
which we call the set of "allowed initial orientations." We combine this set with the best-fit values 

of the remaining measured orbital parameters of planets C and D to define the set of "allowed 
initial conditions." Each member of the set of allowed initial conditions may be used as initial 
conditions for the SWIFT orbit integrator. We refer to the resultant evolution as a "scenario." 

3.2. Sampling 

For each of 19 values of Qq spaced uniformly between 0° and 180°, we randomly generate 
sets of (r2o,XO)0)^) using the following normalized probability distribution functions. Isotropy of 
physical space demands that dP/d<^ = l/27r, dP/d{cosS) = 1/2, and dP/dClo = l/27r; i.e., (p, 
cos 6, and JIq are uniformly distributed over their respective domains. The probability distribution 
function for xoi like that for Oq, is unknown (and may well depend on other parameters such as 
0,q), reflecting the unknown physics of the formation of these planets. We assume for simplicity 
that xo is uniformly distributed over its domain; i.e., dP/dxo = l/27r. 

For a given Qq, those randomly generated values of {0,q,xo,4'i^) which satisfy the observed 
present-day values of uc e [-111?1, -106?9] (la) and ojd e [-106?9, -101?1] (la) [as computed 
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using equation (3)] constitute our sample of allowed initial orientations. For our choice of dP/dxo, 
we find empirically that the fraction of initial orientations that are allowed does not vary with 
Go- At every 0o, the first = 45, randomly generated, allowed initial conditions arc employed 
as input conditions for the SWIFT orbit integrator. We believe that our choice for N adequately 
samples parameter space since we have taken N = 135 at several values of Gq and find none of our 
conclusions below to be changed. As discussed in §2.2, the starting epoch of all integrations was 
JD 2450160.1 days. 



3.3. Dynamics 

Each orbit integration lasts t = 10^ yr. In addition to recording the usual orbital parameters 
such as the eccentricities and mutual inclination, we also compute Ao; according to equation (3). 
The integration spans 10^-10'^ libration/circulation periods of Aw; instabilities usually manifest 
themselves on shorter timescales. A system is considered to be unstable if either planet collides 
with the central star or if either planet is ejected from the system. In practice, we take the 
latter criterion to be satisfied if the semi-major axis of the orbit exceeds 20 AU; we have checked a 
posteriori that stable scenarios preserve semi-major axes to within ~10%. Integrating 45 x 19 = 855 
scenarios, each for t = 10^ yr with a timestep of At = 3 days, requires 6 CPU days on a DEC Alpha 
workstation. 

Secular libration of Ac<j about 0° is a typical outcome for Qq < 30°. Figures 2ghi display the 
evolution of eccentricities, Ao;, and Q in one example scenario for which Qq = 20°. The caption 
to Figure 2 lists the employed values of sin ic and sin i d . The gross features of the evolution are 
the same as those for the co-planar case (cf. Figures 2def). Because the Doppler-measured uj is 
sensitive to both e and 1, Aa; is modulated not just at the single usual frequency but also at 
the two remaining secular precession frequencies, '^g- and the nodal precession frequency 



^ ^ ~ 2M::fe"'V2(")(« + V^) ^ -27r/(6000yr) . (10) 

The above expression is derived to lowest order in Q and the numerical evaluation takes smic ~ 
sin ZD « 1 (as is appropriate for the scenario in Figures 2ghi). In Figure 2i, the mutual inclination 
remains constant at 20°; secular invariance of 6 obtains for 2-planet systems with small and 
small eccentricities (Brouwer & Clemence 1961). 

Circulation of Aw, instability, and, more rarely, libration of Aw are three possible outcomes 
for 40° < 00 90°. The qualitatively distinct character of the disturbing function at such large 
inclinations was first explored by Kozai (1962), who described the secular evolution of highly 
inclined orbits of asteroids in Jupiter's gravitational field. Here planet D, the main repository 
of angular momentum, behaves as Jupiter, and planet C behaves as an asteroid. For mutual 
inclinations such that cos^ < 3/5 (« cos^ 39? 2), the system may inhabit the Kozai resonance, for 



- 12 - 



which Ic and etc precess about the axis parallel to the total angular momentum vector of the system 
at a common frequency while precesses about that axis at the slower frequency ~(7_. The 

result is that Auj circulates at ^\^\- Figures 2jkl show how the eccentricities, Auj, and Q evolve 
under the Kozai mechanism. Note how ec and G undergo larger variations than in the small Q case; 
ec and 6 can be coupled in the cos^ G < 3/5 case such that both vary secularly at frequency ~2|r2| 
while keeping Kozai's integral, ^1 — e^cosQ, approximately constant. The secular driving of ec 
to values approaching unity causes 90% of our sampled scenarios having Go ~ 90° to self-destruct 
within t = 10^ yr. 

For Go > 90°, the angular momenta of planets C and D, projected onto the total angular 
momentum vector of the system, point in opposite directions. Consequently, secular exchanges of 
angular momenta between the two planets cause their eccentricities to rise and fall together (Rivera 
& Lissauer 2000a). For 90° ^ Gq ^ 140°, the synchronized oscillations of the eccentricities com- 
pounds the destabilizing effects of coupled ec and G to render the majority of scenarios unstable. 
For 150° < Go < 180°, ec and G decouple and Aa; circulates stably. A typical example of the 
evolution in the nearly co-planar, counter-rotating case is showcased in Figures 2mno. 

Figure 4 summarizes the effects we have discussed so far. The solid line traces the fraction 
of allowed initial conditions, sampled according to the assumptions stated in §3.2, that generate 
2-planet scenarios which survive for t = 10^ yr. Nearly co-planar scenarios, whether co-rotating 
or counter-rotating, betray little or no sign of instability. Those co-rotating systems at Gq < 30° 
that are unstable all have sini < 0.3 for either planet C or D. Our estimated fraction of surviving 
systems is an upper limit because we have neglected thus far the presence of planet B and because 
there may be instabilities that require t > 10^ yr to develop. We expect the survival fraction to be 
most severely overestimated at 40° < Go < 140°, inclinations for which ec and G undergo large, 
coupled variations; values of ec > 0.93 can cause planets C and B to undergo destabilizing close 
encounters. 

Open circles in Figure 4 mark Pume (l^f^l < 10°), the fraction of time a surviving sampled 
scenario spends with \Auj\ < 10°. For Gq ^ 20°, the two planets are locked in a secular apsidal 
resonance for which Pume can be as high as 0.3-0.4. For Gq ^ 40°, Alu typically circulates so that 
Ptime = 20°/360° = 0.055. The specific distribution of values of Pume at a given Gq depends upon 
our assumed sampling function for xoj the initial argument of pericenter (see §3.2). For given Gq, 
lower values of Pume typically correspond to lower values of sinic and sinin (see §2.2). Maximum 
values of Pume ~ 0.39 are attained at Go = 10-20°. 

Surprisingly, a secular apsidal resonance exists at Go ~ 70-80° for which Pume ~ 0.2. The 
evolution within this "anomalous" resonance at high Go is characterized by rapid oscillations in 
Au, ec, eD, and G at twice the nodal precession frequency, ~2|i7|. The eccentricity of planet C 
is periodically driven by secular perturbations to values as high as 0.95. Given such large values 
of ec, we suspect that including the presence of planet B may render these anomalously resonant 
configurations unstable. We confirm this suspicion in the next section. 
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4. Planets B, C, and D 

In the previous section we found those allowed initial conditions involving only planets C and 
D that result in libration of Aw. Here we add the effects of the third, innermost planet B to these 
resonant configurations, in the expectation that some of the resultant 3-planet systems will not be 
stable. 

The inclusion of planet B introduces a number of complications. First, three extra dimensions 
are added to the space of allowed initial orientations: Qbo, ^bo, and xbo — respectively, the initial 
inclination, initial longitude of ascending node, and initial argument of pericenter of the orbit of 
planet B upon the reference orbit of planet C (see §3.1). Since the present-day orbit of planet 
B appears nearly indistinguishable from a circle (see Table 1), we can take e^o = to eliminate 
Xbo as a dimension. Even with this simplifying assumption, the volume of additional parameter 
space to be surveyed is, in principle, substantial. Computational needs are further exacerbated by 
the fact that the orbital period of planet B is 4.6 days; this necessitates a commensurately short 
computational timestep. Finally, planet B is sufficiently close to its parent star that contributions 
to the planet's apsidal and nodal precession rates due to stellar oblateness and general relativistic 
effects are not negligible compared to contributions from the outer two planets. 

We proceed with a simplified program. First, we select a set of allowed initial orientations for 
planets C and D which give rise to apsidally resonant, 2-planet scenarios. At each Gq, we select 
the one 2-planet scenario that generates the maximum value of Pume d^"^! < 10°), provided this 
quantity exceeds 0.15. As evident in Figure 4, the values of 9o which give rise to such scenarios are 
0°, 10°, 20°, 30°, 40°, 70°, and 80°. Next, to each member of this set of allowed initial orientations, 
we add the following orientation angles for planet B: Qbo = ©o and Qbo = ^o- That is, the 
orbit of planet B is assumed to be initially circular and to lie in the orbit plane of the most massive 
planet, D. Finally, this augmented set of allowed initial orientations is combined with the remaining 
measured orbital parameters for all three planets and employed as initial conditions for the SWIFT 
orbit integrator. We employ a timestep of At = 0.20 days (5% of the orbital period of planet B) 
and run each integration for i = 10^ yr. The effects of stellar oblateness and general relativity are 
ignored. 

While the above procedure does not permit accurate determination of the long-term orbital 

evolution of planet B, we believe it does capture qualitatively the effect of planet B on the possibility 
of resonant apsidal lock between planets C and D. We expect the innermost planet to represent 
primarily a potential means of destroying the apsidal resonance, as the secular growth of ec takes 
planet C disruptively close to planet B. 

Scenarios for which all three planets survive the entire duration of the integration are shown 
in Figure 5. Only scenarios for which 6o = ©so < 40° qualify. For 6o > 70°, the eccentricity 
of planet C is secularly driven to such high values that planet B is perturbed by planet C into a 
star-crossing orbit in f ~ 10^ yr. Thus, the "anomalous" resonant configurations at @o = 70°-80° 
found in §3.3 do not survive the inclusion of planet B. A similar fate probably befalls the scenario 
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for which 9o = 40° ; an orbit integration of duration t > 10^ yr is required to confirm this suspicion. 

For 00 < 30° , planets C and D can remain in stable resonant lock despite the presence of planet 
B. Values of Pume (l^^l < 10°) in 3-planet scenarios are somewhat reduced from their respective 
values in 2-planet scenarios. Values of Pume ~ 0.28-0.33 obtain for 6o < 20°; at 0o = 30°, Pume 
drops to 0.13. 

5. Concluding Remarks 

Our results suggest that the outer two planets of Upsilon Andromedae occupy nearly co-planar, 
co-rotating orbits. The lines of evidence are as follows: 

1. Finding the Doppler measured w's of planets C and D to be nearly equal today is least 
surprising if the mutual orbit inclination 0o < 20° (Figures 4 and 5) . If 9o < 20° , planets C 
and D can be locked in a secular apsidal resonance for which Au = ojd — ojc librates about 0° 
with an amplitude of approximately 30°. The corresponding fraction of time that | Awl < 10° 
is Ptime ^ 0.30. These libration parameters are subject to change with future refinements 
in the values of the present-day fitted eccentricities (Figure 3). By contrast, if Oq ^ 40°, 
typically Aw circulates and Pume ~ 20°/360° = 0.055 for the relatively few allowed initial 
conditions that give rise to stable systems. 

2. The secular apsidal resonance at Qq < 20° limits variations in eccentricities and thereby 
affords the system stability (Figure 2). By contrast, if 40° < 0o ^ 90°, secular driving 
of the eccentricity of planet C to values approaching unity via the Kozai resonance fosters 
close approaches and concomitant instability. If 90° < 0o 140°, the destabilizing effects of 
coupled ec and are abetted by synchronized driving of the eccentricities. 

We note that co-planar, counter-rotating configurations (0 ~ 180°) are allowed by argument 2 but 
disfavored by argument 1. Furthermore, values of sinic, sini^i > 0.5 are favored, not only for 
reasons of stability, but also because such values are typical of scenarios for which the fraction of 
time spent with small |Aa;| is maximized. 

The preference for small mutual inclination, based on combining considerations of stability 
with the likelihood of finding |Aa;| <C Irad, is consistent with the origin of these planets in a 
flattened, circumstellar disk. How does our expectation that ©o < 20° compare with orbital 

inclinations in related contexts? The opening angle of hydrostatically flared protoplanetary disks 
at stellocentric distances of '^l AU is ~6° (see, e.g., Chiang & Goldreich 1997). The planets in the 
Solar System, with the understandable exception of Pluto (Malhotra 1998), occupy orbits whose 
mutual inclinations do not exceed about 7 degrees; Mercury aside, mutual inclinations do not 
exceed about 3 degrees. 
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By contrast, scenarios for the formation of extrasolar planets that do not involve a primordial 
disk may predict a substantial population of systems for which > 150°. For example, Papaloizou 
& Terquem (2000) have proposed that extrasolar planets might gravitationally fragment from, and 
scatter out of, a circumstellar spherical cloud of gas. Planetary systems having a wide distribution 
of initial mutual inclinations from @o = 0° to 180° might be expected to form. Prom the stability 
considerations elucidated above, we might expect a bimodal population of systems to survive — 
those for which 0° < < 40° and those for which 150° < B < 180° — intermediate inclinations 
having been eliminated via the destructive Kozai effect. 

The final word on the actual values of orbital inclinations in extrasolar planetary systems will 
likely come from combining the Doppler radial velocity data with astrometric measurements, or, 
in rare cases, from transit observations (Mazeh et al. 2000; Queloz et al. 2000). The amplitudes 
of stellar wobbles induced by planets C and D are (0.0943/ sinic) mas and (0.576/ sini/)) mas, 
respectively (Pourbaix 2001). The star v And is sufficiently bright {V = 4.2) that the upcoming 
astrometric satellite FAME (Full-Sky Astrometric Mapping Explorer) should yield at least a 2a 
detection of the astrometric signature of planet C. 

We thank Geoff Marcy and Debra Fischer for generously providing updated values of the 
fitted orbital parameters of v And, Fathi Namouni for useful conversations, and Greg Laughlin 
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Fig. 1. — Defining geometry for the Doppler-measured argument of pericenter u, the angle between 
pericenter and the node of the planetary orbit on the plane of the sky. The sky plane selected is 
that which passes through the system barycenter. The vector n lies in the planes of the orbit and 
of the sky. The vector e lies in the orbit plane. The node selected is the one for which the planet 
travels towards the observer. The argument u advances along the orbit plane in the direction of 
increasing true anomaly. 

Fig. 2. — Sampling of 2-planet scenarios in which the orbits of planets C and D are initially 
inclined by the value of 9o indicated. Each horizontal row of panels displays the evolution of ec, 
e£), Aw = ojD — ijJc^ and = arccos (Ic'Id) for a given scenario. For panels (a)-(c), the evolution is 
computed using the classical, analytic solution of Laplace-Lagrangc for B = 0° . For the remaining 
panels, the evolution is computed by numerical integration. Libration of AtJ is a typical outcome 
for 00 ^ 20°, while circulation of Acj and large secular variations of ec and Q via the Kozai 
mechanism are possible outcomes for 40° < Oo ^ 140°. For Oq > 150°, Ao; circulates while ec 
and vary modestly and synchronously. Note in panel (f) how Aa; is modulated at a number 
of frequencies, refiecting its dependence on the 4 precessing vectors (l)c,D and (e)c,D- For panels 
(a)-(f), sini = 1 and timcscalcs scale as sini. For panels (g)-(i), sinic = 0.99 and sini/j = 0.93. 
For panels (j)-(l), sinic = 1 and smio = 0.73. For panels (m)-(o), sinic = 0.78 and smir) = 0.87. 

Fig. 3. — Contours of constant libration amplitude in degrees as a function of possible present-day 
fitted eccentricities, for = 0° and sini = 1. The cross indicates the best-fit eccentricities from 
Marcy &: Fischer (2001), while the inner and outer dotted ellipses enclose the la and 2a confidence 
regions for the eccentricities, respectively. Panel (a) is calculated using Laplace-Lagrange theory 
while panel (b) plots the maximum value of Aa; over t = 10^ yr in numerical integrations using 
SWIFT. 

Fig. 4. — Open circles, left-hand ordinate: Fraction of time that surviving, sampled, 2-planet 
scenarios spend with |Au;| < 10°. The observation today that |Au;| < 10° is least surprising if 
the mutual inclination between the orbit planes of planets C and D is less than or equal to 
20°. At these modest inclinations, the two planets inhabit a secular apsidal resonance in which Ao; 
librates about 0° with an amplitude of ~25°. The cluster of points at = 70°-80° and relatively 
large Pume ~ 0.20 refiects the existence of an "anomalous" apsidal 2-planet resonance. In §4, we 
discuss how the inclusion of the third, innermost planet B disrupts the anomalous resonance. Solid 
line, right-hand ordinate: Fraction of 2-planet scenarios that survive the 10^ yr duration of the 
integration. Nearly co-planar orbits, whether co-rotating or counter-rotating, betray little sign of 
instability. Scenarios at intermediate inclinations tend to be unstable because the Kozai effect at 
> 40°, abetted by synchronous driving of eccentricities at > 90°, drives the eccentricity of 
planet C to such high values that planet C collides with the star. 

Fig. 5. — Three-planet scenarios in which planets B and D initially occupy co-planar orbits while 
planet C initially occupies an orbit inclined by the value of 0o indicated. Each horizontal row 
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of panels displays the evolution of ec, eo, Aa; = ud — ^Ci and = arccos (Ic - Id) for a given 
scenario. Scenarios for which ©q < 20° preserve apsidal libration of Au; despite the presence 
of planet B and spend the maximum fraction of time with |Aa;| < 10°: Pume ^ 0.28-0.33. At 
Oo = 30°, the amplitude of apsidal libration is at its largest, approximately 90°. At Oq > 40°, 
Aa; eventually circulates. For panels (a)-(c), sini^ = sinic = sinij) = 0.89. For panels (d)- 
(f), sinzfi = sinzD = 0.99 and sinzc = 0.95. For panels (g)-(i), sinie = sinzo = 0.79 and 
smic = 0.55. For panels (j)-(l), sini^ = sini/j = 0.99 and sinic = 0.97. For panels (m)-(o), 
siiLiB = siiLiD = 0.90 and sinic = 0.44. 



Table 1. Fitted Orbital Parameters of Upsilon Andromedae (Marcy & Fischer 2001) 



Planet 


P (days) 




e 


a; 


K (m/s) 


msinz (Mj) 


a (AU) 


B 


4.61706 (7 X 10-^) 


2450001.6 (NA)'= 


0.015 (0.009) 


32.1 (NA)^ 


71.0 (0.7) 


0.69 


0.059 


C 


241.14 (0.22) 


2450160.1 (1.9) 


0.254 (0.015) 


-109.0 (2.9) 


55.6 (0.9) 


1.96 


0.828 


D 


1309.13 (5.13) 


2450044.0 (11.3) 


0.342 (0.015) 


-104.2 (3.8) 


66.2 (1.2) 


3.98 


2.557 



^Based on fitting A'" = 189 radial velocity points, with rms residuals of 11.5 m/s. 

''The orbital period is P, the Julian date of pericenter passage is TpeH, the eccentricity is e, the planetary mass 
is m measured in Jupiter masses Mj, and the semi-major axis is a. Remaining variables are defined in the text. 
Uncertainties (lo") in fitted values arc enclosed in parentheses. 

■^Uncertainties in u) and Tperi have little meaning for the nearly circular orbit of planet B. 



